7. Plotting primary analyses
In these plots we summarize the 100 partial dependence plots run for each analysis as envelopes. The lighter gray envelope describes the maximum and minimum values for each timestep and the darker gray envelope describes the 25th and 75th percentiles.
Molybdenum
mo.pdp <- interp.pdp(data = Mo.anox.py.rf.partial, age.rounding.factor = 1, n = n)
Mo.plot <- ggplot()+
geom_ribbon(data = mo.pdp, aes(x = Age, ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(data = mo.pdp, aes(x = Age, ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(0,62),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Mo (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.3,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
Mo.plot

Uranium
u.pdp <- interp.pdp(data = U.anox.rf.partial, age.rounding.factor = 1, n = n)
U.plot <- ggplot()+
geom_ribbon(data = u.pdp, aes(x = Age, ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(data = u.pdp, aes(x = Age, ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(0,62),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("U (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.3,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
U.plot

Proportion euxinic
Note that we hide warnings here because of random forest warning about using binary variable in a regression (here we do want to do a regression to replicate the broad approach of Sperling et al. 2015, Nature - see similar comments above).
prop_eux.pdp <- interp.pdp(data = Fepy.anox.rf.partial, age.rounding.factor = 1, n = n)
prop_eux.plot <- ggplot()+
geom_ribbon(data = prop_eux.pdp, aes(x = Age, ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(data = prop_eux.pdp, aes(x = Age, ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-.01,1.03),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Proportion euxinic")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.3,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
prop_eux.plot

Total Organic Carbon
TOC.pdp <- interp.pdp(data = TOC.all.rf.partial, age.rounding.factor = 1, n = n)
TOC.plot <- ggplot()+
geom_ribbon(data = TOC.pdp, aes(x = Age, ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(data = TOC.pdp, aes(x = Age, ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(0,8.3),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("TOC (wt %)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.3,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
TOC.plot

Combine to generate primary summary plot
In these plots we have delineated the three main phases of marine biogeochemistry between 1000 and 300Ma - with transitions around the base of of the Cambrian and Devonian periods.
Mo.plot.for.sum <- ggplot(mo.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.4,62),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Mo (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
U.plot.for.sum <- ggplot(u.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-.4,62),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("U (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
prop_eux.plot.for.sum <- ggplot(prop_eux.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.01,1.03),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Proportion euxinic")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
TOC.plot.for.sum <- ggplot(TOC.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-0.04,8.3),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("TOC (wt%)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
summary.plot <- ggarrange2(Mo.plot.for.sum, prop_eux.plot.for.sum, U.plot.for.sum, TOC.plot.for.sum, ncol=2, heights=c(1,1))

ggsave("Figure 2 Partial Dependence Plot with shading.pdf", summary.plot, height=14, width=26)
8. Expanded analyses varying redox-sensitive predictors
Anoxic Mo (no pyrite)
Here, we generate an anoxic Mo equivalent to our primary U dataset (not restricting to anoxic samples with iron speciation). Fepy/FeHR is not included as a predictor variable.
Mo.anox.rf <- trace.toc.full %>%
filter(!is.na(Mo..ppm.)) %>%
filter(FeHR.FeT >= 0.38 | FeT.Al >= 0.53)
nrow(Mo.anox.rf)
Mo.anox.rf.results <- Monte.Carlo.rF(data = Mo.anox.rf,
resp.var = "Mo..ppm.",
vars = c("age.model",
"site.type",
"metamorphic.bin",
"basin.type",
"site.latitude",
"site.longitude",
"lithology.name",
"environmental.bin",
"Mo..ppm.",
"TOC..wt..",
"Al..wt.."),
n = n,
run.age.unc = TRUE,
run.partial.context = TRUE,
run.Al.impute = TRUE)
Mo.anox.rf.partial <- Mo.anox.rf.results$partial.plot.data
Mo.anox.rf.import <- Mo.anox.rf.results$importance.data
mo.no_py.pdp <- interp.pdp(data = Mo.anox.rf.partial, age.rounding.factor = 1, n = n)
We combine these analyses with our primary Mo analyses for plotting.
mo.pdp$treatment <- "Anoxic samples + TOC + Fepy/FeHR"
mo.no_py.pdp$treatment <- "Anoxic samples + TOC"
mo.pdp.sum <- rbind(mo.pdp,
mo.no_py.pdp)
Anoxic U with pyrite
Here, we generate an anoxic U equivalent to our primary Mo dataset (restricting to anoxic samples with iron speciation). Fepy/FeHR is included as a predictor variable.
U.anox.py.rf <- trace.toc.full %>%
filter(!is.na(U..ppm.) & !is.na(Fe.py.FeHR)) %>%
filter(FeHR.FeT >= 0.38)
nrow(U.anox.py.rf)
U.anox.py.rf.results <- Monte.Carlo.rF(data = U.anox.py.rf,
resp.var = "U..ppm.",
vars = c("age.model",
"site.type",
"metamorphic.bin",
"basin.type",
"site.latitude",
"site.longitude",
"lithology.name",
"environmental.bin",
"U..ppm.",
"TOC..wt..",
"Fe.py.FeHR",
"Al..wt.."),
n = n,
run.age.unc = TRUE,
run.partial.context = TRUE,
run.Al.impute = TRUE)
U.anox.py.rf.partial <- U.anox.py.rf.results$partial.plot.data
U.anox.py.rf.import <- U.anox.py.rf.results$importance.data
u.w_py.pdp <- interp.pdp(data = U.anox.py.rf.partial, age.rounding.factor = 1, n = n)
We combine these analyses with our primary U analyses for plotting.
u.pdp$treatment <- "Anoxic samples + TOC"
u.w_py.pdp$treatment <- "Anoxic samples + TOC + Fepy/FeHR"
u.pdp.sum <- rbind(u.pdp,
u.w_py.pdp)
Proportion euxinic with TOC
Note that we hide warnings here because of random forest warning about using binary variable in a regression (here we do want to do a regression to replicate the broad approach of Sperling et al. 2015, Nature).
Fepy.anox.w_TOC.rf <- trace.toc.full %>%
filter(!is.na(Fe.py.FeHR) & !is.na(TOC..wt..)) %>%
filter(FeHR.FeT >= 0.38)
nrow(Fepy.anox.w_TOC.rf)
# For the analysis of the proportion of euxinic samples, we also need to code samples based # upon whether they are euxinic (based on iron speciation) in a binary fashion.
Fepy.anox.w_TOC.rf$euxinic.Fe[Fepy.anox.w_TOC.rf$Fe.py.FeHR >= 0.7] <- 1
Fepy.anox.w_TOC.rf$euxinic.Fe[Fepy.anox.w_TOC.rf$Fe.py.FeHR < 0.7] <- 0
Fepy.anox.w_TOC.rf.results <- Monte.Carlo.rF(data = Fepy.anox.w_TOC.rf,
resp.var = "euxinic.Fe",
vars = c("age.model",
"site.type",
"metamorphic.bin",
"basin.type",
"site.latitude",
"site.longitude",
"lithology.name",
"environmental.bin",
"euxinic.Fe",
"Al..wt..",
"TOC..wt.."),
n = n,
run.age.unc = TRUE,
run.partial.context = TRUE,
run.Al.impute = TRUE)
Fepy.anox.w_TOC.rf.partial <- Fepy.anox.w_TOC.rf.results$partial.plot.data
Fepy.anox.w_TOC.rf.import <- Fepy.anox.w_TOC.rf.results$importance.data
prop_eux.w_TOC.pdp <- interp.pdp(data = Fepy.anox.w_TOC.rf.partial, age.rounding.factor = 1, n = n)
We combine these analyses with our primary proportion euxinic analyses for plotting.
prop_eux.pdp$treatment <- "Anoxic samples"
prop_eux.w_TOC.pdp$treatment <- "Anoxic samples + TOC"
prop_eux.pdp.sum <- rbind(prop_eux.pdp,
prop_eux.w_TOC.pdp)
Anoxic TOC (no pyrite)
TOC.anox.rf <- trace.toc.full %>%
filter(!is.na(TOC..wt..)) %>%
filter(FeHR.FeT >= 0.38 | FeT.Al >= 0.53)
nrow(TOC.anox.rf)
TOC.anox.rf.results <- Monte.Carlo.rF(data = TOC.anox.rf,
resp.var = "TOC..wt..",
vars = c("age.model",
"site.type",
"metamorphic.bin",
"basin.type",
"site.latitude",
"site.longitude",
"lithology.name",
"environmental.bin",
"TOC..wt..",
"Al..wt.."),
n = n,
run.age.unc = TRUE,
run.partial.context = TRUE,
run.Al.impute = TRUE)
TOC.anox.rf.partial <- TOC.anox.rf.results$partial.plot.data
TOC.anox.rf.import <- TOC.anox.rf.results$importance.data
TOC.anox.pdp <- interp.pdp(data = TOC.anox.rf.partial, age.rounding.factor = 1, n = n)
Anoxic TOC with pyrite
TOC.anox.py.rf <- trace.toc.full %>%
filter(!is.na(TOC..wt..) & !is.na(Fe.py.FeHR)) %>%
filter(FeHR.FeT >= 0.38)
nrow(TOC.anox.rf)
TOC.anox.py.rf.results <- Monte.Carlo.rF(data = TOC.anox.py.rf,
resp.var = "TOC..wt..",
vars = c("age.model",
"site.type",
"metamorphic.bin",
"basin.type",
"site.latitude",
"site.longitude",
"lithology.name",
"environmental.bin",
"TOC..wt..",
"Al..wt..",
"Fe.py.FeHR"),
n = n,
run.age.unc = TRUE,
run.partial.context = TRUE,
run.Al.impute = TRUE)
TOC.anox.py.rf.partial <- TOC.anox.py.rf.results$partial.plot.data
TOC.anox.py.rf.import <- TOC.anox.py.rf.results$importance.data
TOC.anox.py.pdp <- interp.pdp(data = TOC.anox.py.rf.partial, age.rounding.factor = 1, n = n)
We combine these analyses with our primary TOC analyses for plotting.
TOC.pdp$treatment <- "All samples"
TOC.anox.pdp$treatment <- "Anoxic samples"
TOC.anox.py.pdp$treatment <- "Anoxic samples + Fepy/FeHR"
TOC.pdp.sum <- rbind(TOC.pdp,
TOC.anox.pdp,
TOC.anox.py.pdp)
Summary plotting of analyses including varying redox treatments
In these summary plots we just plot the interquartile ranges for our analyses to ease the comparison of trends in central tendancy.
Mo.plot.for.redox.var.sum <- ggplot(mo.pdp.sum, aes(x=Age))+
#annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75, fill=treatment, color=treatment), alpha=.7, size=.5)+
scale_fill_manual(values = c(rgb(127, 160, 86, maxColorValue = 255),
rgb(179, 224, 149, maxColorValue = 255)),
name = "Treatment")+
scale_color_manual(values = c(rgb(127, 160, 86, maxColorValue = 255),
rgb(179, 224, 149, maxColorValue = 255)),
name = "Treatment")+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.4,82),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,254,600,700,800,900,1000))+
ylab("Mo (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
U.plot.for.redox.var.sum <- ggplot(u.pdp.sum, aes(x=Age))+
#annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75, fill=treatment, color=treatment), alpha=.7, size=.5)+
scale_fill_manual(values = c(rgb(103, 165, 153, maxColorValue = 255),
rgb(191, 208, 186, maxColorValue = 255)),
name = "Treatment")+
scale_color_manual(values = c(rgb(103, 165, 153, maxColorValue = 255),
rgb(191, 208, 186, maxColorValue = 255)),
name = "Treatment")+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-.4,62),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("U (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
prop_eux.plot.for.redox.var.sum <- ggplot(prop_eux.pdp.sum, aes(x=Age))+
#annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75, fill=treatment, color=treatment), alpha=.7, size=.5)+
scale_fill_manual(values = c(rgb(240, 64, 40, maxColorValue = 255),
rgb(251, 154, 133, maxColorValue = 255)),
name = "Treatment")+
scale_color_manual(values = c(rgb(240, 64, 40, maxColorValue = 255),
rgb(251, 154, 133, maxColorValue = 255)),
name = "Treatment")+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.01,1.03),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Proportion euxinic")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
TOC.plot.for.redox.var.sum <- ggplot(TOC.pdp.sum, aes(x=Age))+
#annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75, fill=treatment, color=treatment), alpha=.7, size=.5)+
scale_fill_manual(values = c(rgb(203, 140, 55, maxColorValue = 255),
rgb(254, 179, 66, maxColorValue = 255),
rgb(254, 217, 106, maxColorValue = 255)),
name = "Treatment")+
scale_color_manual(values = c(rgb(203, 140, 55, maxColorValue = 255),
rgb(254, 179, 66, maxColorValue = 255),
rgb(254, 217, 106, maxColorValue = 255)),
name = "Treatment")+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-0.04,8.3),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("TOC (wt%)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
summary.redox.var.plot <- ggarrange2(Mo.plot.for.redox.var.sum, prop_eux.plot.for.redox.var.sum, U.plot.for.redox.var.sum, TOC.plot.for.redox.var.sum, ncol=2, heights=c(1,1))

ggsave("Figure Sx Partial Dependence Plot (varying redox treatments).pdf", summary.redox.var.plot, height=14, width=26)
9. Primary analyses without Al
We also re-run our analyses without incorporating [Al] as a broad proxy for detrital input.
Molybdenum
Mo.anox.py.no_Al.rf.results <- Monte.Carlo.rF(data = Mo.anox.py.rf,
resp.var = "Mo..ppm.",
vars = c("age.model",
"site.type",
"metamorphic.bin",
"basin.type",
"site.latitude",
"site.longitude",
"lithology.name",
"environmental.bin",
"Mo..ppm.",
"TOC..wt..",
"Fe.py.FeHR"),
n = n,
run.age.unc = TRUE,
run.partial.context = TRUE,
run.Al.impute = FALSE)
Mo.anox.py.no_Al.rf.partial <- Mo.anox.py.no_Al.rf.results$partial.plot.data
Mo.anox.py.no_Al.rf.import <- Mo.anox.py.no_Al.rf.results$importance.data
mo.no_Al.pdp <- interp.pdp(data = Mo.anox.py.no_Al.rf.partial, age.rounding.factor = 1, n = n)
Uranium
U.anox.no_Al.rf.results <- Monte.Carlo.rF(data = U.anox.rf,
resp.var = "U..ppm.",
vars = c("age.model",
"site.type",
"metamorphic.bin",
"basin.type",
"site.latitude",
"site.longitude",
"lithology.name",
"environmental.bin",
"U..ppm.",
"TOC..wt.."),
n = n,
run.age.unc = TRUE,
run.partial.context = TRUE,
run.Al.impute = FALSE)
U.anox.no_Al.rf.partial <- U.anox.no_Al.rf.results$partial.plot.data
U.anox.no_Al.rf.import <- U.anox.no_Al.rf.results$importance.data
u.no_Al.pdp <- interp.pdp(data = U.anox.no_Al.rf.partial, age.rounding.factor = 1, n = n)
Proportion euxinic
Note that we hide warnings here because of random forest warning about using binary variable in a regression (here we do want to do a regression to replicate the broad approach of Sperling et al. 2015, Nature).
Fepy.anox.no_Al.rf.results <- Monte.Carlo.rF(data = Fepy.anox.rf,
resp.var = "euxinic.Fe",
vars = c("age.model",
"site.type",
"metamorphic.bin",
"basin.type",
"site.latitude",
"site.longitude",
"lithology.name",
"environmental.bin",
"euxinic.Fe"),
n = n,
run.age.unc = TRUE,
run.partial.context = TRUE,
run.Al.impute = FALSE)
Fepy.anox.no_Al.rf.partial <- Fepy.anox.no_Al.rf.results$partial.plot.data
Fepy.anox.no_Al.rf.import <- Fepy.anox.no_Al.rf.results$importance.data
prop_eux.no_Al.pdp <- interp.pdp(data = Fepy.anox.no_Al.rf.partial, age.rounding.factor = 1, n = n)
Total Organic Carbon
TOC.all.no_Al.rf.results <- Monte.Carlo.rF(data = TOC.all.rf,
resp.var = "TOC..wt..",
vars = c("age.model",
"site.type",
"metamorphic.bin",
"basin.type",
"site.latitude",
"site.longitude",
"lithology.name",
"environmental.bin",
"TOC..wt.."),
n = n,
run.age.unc = TRUE,
run.partial.context = TRUE,
run.Al.impute = FALSE)
TOC.all.no_Al.rf.partial <- TOC.all.no_Al.rf.results$partial.plot.data
TOC.all.no_Al.rf.import <- TOC.all.no_Al.rf.results$importance.data
TOC.no_Al.pdp <- interp.pdp(data = TOC.all.no_Al.rf.partial, age.rounding.factor = 1, n = n)
Combine to generate summary plot
Mo.no_Al.plot.for.sum <- ggplot(mo.no_Al.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.4,62),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Mo (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
U.no_Al.plot.for.sum <- ggplot(u.no_Al.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-.4,62),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("U (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
prop_eux.no_Al.plot.for.sum <- ggplot(prop_eux.no_Al.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.01,1.03),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Proportion euxinic")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
TOC.no_Al.plot.for.sum <- ggplot(TOC.no_Al.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-0.04,8.3),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("TOC (wt%)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
summary.plot <- ggarrange2(Mo.no_Al.plot.for.sum, prop_eux.no_Al.plot.for.sum, U.no_Al.plot.for.sum, TOC.no_Al.plot.for.sum, ncol=2, heights=c(1,1))

ggsave("Figure Sx Partial Dependence Plot (no Al) with shading.pdf", summary.plot, height=14, width=26)
10. Variable importance plots
We will summarize the variable importance of the predictor variables in our primary Monte Carlo random forest analyses, using box and whisker plots to illustrate the distributions generated by our Monte Carlo approach.
Mo.MSE.plot <- ggplot(Mo.anox.py.rf.import, aes(x = Variable, y = IncMSE))+stat_boxplot(geom ='errorbar', width = 0.6, size=0.6, color="grey20", lwd=.1)+ geom_boxplot(outlier.alpha = .5, outlier.shape=16, outlier.size=2, outlier.color="grey20", color="grey20", lwd=0.6, fill="grey80", fatten = 1.3)+coord_flip()+
theme_bw()+
ylab("Increased Mean\nSquared Error (%)")+
ggtitle("Molybdenum")+
theme(panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1.1),
axis.title = element_text(size=30),
axis.text = element_text(size=24, colour="black"),
plot.title = element_text(size=30, face = "bold"),
plot.margin = ggplot2::margin(10,30,10,10),
legend.position="none",
axis.title.y = element_blank(),
axis.title.x = element_text(margin = ggplot2::margin(t = 10, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
## Warning: Duplicated aesthetics after name standardisation: size
Mo.Node.plot <- ggplot(Mo.anox.py.rf.import, aes(x = Variable, y = IncNodePurity))+stat_boxplot(geom ='errorbar', width = 0.6, size=0.6, color="grey20", lwd=.1)+ geom_boxplot(outlier.alpha = .5, outlier.shape=16, outlier.size=2, outlier.color="grey20", color="grey20", lwd=0.6, fill="grey80", fatten = 1.3)+coord_flip()+
theme_bw()+
ylab("Increased\nNode Purity")+
ggtitle("Molybdenum")+
theme(panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1.1),
axis.title = element_text(size=30),
axis.text = element_text(size=24, colour="black"),
plot.title = element_text(size=30, face = "bold"),
plot.margin = ggplot2::margin(10,30,10,10),
legend.position="none",
axis.title.y = element_blank(),
axis.title.x = element_text(margin = ggplot2::margin(t = 10, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
## Warning: Duplicated aesthetics after name standardisation: size
U.MSE.plot <- ggplot(U.anox.rf.import, aes(x = Variable, y = IncMSE))+stat_boxplot(geom ='errorbar', width = 0.6, size=0.6, color="grey20", lwd=.1)+ geom_boxplot(outlier.alpha = .5, outlier.shape=16, outlier.size=2, outlier.color="grey20", color="grey20", lwd=0.6, fill="grey80", fatten = 1.3)+coord_flip()+
theme_bw()+
ylab("Increased Mean\nSquared Error (%)")+
ggtitle("Uranium")+
theme(panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1.1),
axis.title = element_text(size=30),
axis.text = element_text(size=24, colour="black"),
plot.title = element_text(size=30, face = "bold"),
plot.margin = ggplot2::margin(10,30,10,10),
legend.position="none",
axis.title.y = element_blank(),
axis.title.x = element_text(margin = ggplot2::margin(t = 10, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
## Warning: Duplicated aesthetics after name standardisation: size
U.Node.plot <- ggplot(U.anox.rf.import, aes(x = Variable, y = IncNodePurity))+stat_boxplot(geom ='errorbar', width = 0.6, size=0.6, color="grey20", lwd=.1)+ geom_boxplot(outlier.alpha = .5, outlier.shape=16, outlier.size=2, outlier.color="grey20", color="grey20", lwd=0.6, fill="grey80", fatten = 1.3)+coord_flip()+
theme_bw()+
ylab("Increased\nNode Purity")+
ggtitle("Uranium")+
theme(panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1.1),
axis.title = element_text(size=30),
axis.text = element_text(size=24, colour="black"),
plot.title = element_text(size=30, face = "bold"),
plot.margin = ggplot2::margin(10,30,10,10),
legend.position="none",
axis.title.y = element_blank(),
axis.title.x = element_text(margin = ggplot2::margin(t = 10, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
## Warning: Duplicated aesthetics after name standardisation: size
prop_eux.MSE.plot <- ggplot(Fepy.anox.rf.import, aes(x = Variable, y = IncMSE))+stat_boxplot(geom ='errorbar', width = 0.6, size=0.6, color="grey20", lwd=.1)+ geom_boxplot(outlier.alpha = .5, outlier.shape=16, outlier.size=2, outlier.color="grey20", color="grey20", lwd=0.6, fill="grey80", fatten = 1.3)+coord_flip()+
theme_bw()+
ylab("Increased Mean\nSquared Error (%)")+
ggtitle("Proportion Euxinic")+
theme(panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1.1),
axis.title = element_text(size=30),
axis.text = element_text(size=24, colour="black"),
plot.title = element_text(size=30, face = "bold"),
plot.margin = ggplot2::margin(10,30,10,10),
legend.position="none",
axis.title.y = element_blank(),
axis.title.x = element_text(margin = ggplot2::margin(t = 10, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
## Warning: Duplicated aesthetics after name standardisation: size
prop_eux.Node.plot <- ggplot(Fepy.anox.rf.import, aes(x = Variable, y = IncNodePurity))+stat_boxplot(geom ='errorbar', width = 0.6, size=0.6, color="grey20", lwd=.1)+ geom_boxplot(outlier.alpha = .5, outlier.shape=16, outlier.size=2, outlier.color="grey20", color="grey20", lwd=0.6, fill="grey80", fatten = 1.3)+coord_flip()+
theme_bw()+
ylab("Increased\nNode Purity")+
ggtitle("Proportion Euxinic")+
theme(panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1.1),
axis.title = element_text(size=30),
axis.text = element_text(size=24, colour="black"),
plot.title = element_text(size=30, face = "bold"),
plot.margin = ggplot2::margin(10,30,10,10),
legend.position="none",
axis.title.y = element_blank(),
axis.title.x = element_text(margin = ggplot2::margin(t = 10, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
## Warning: Duplicated aesthetics after name standardisation: size
TOC.MSE.plot <- ggplot(TOC.all.rf.import, aes(x = Variable, y = IncMSE))+stat_boxplot(geom ='errorbar', width = 0.6, size=0.6, color="grey20", lwd=.1)+ geom_boxplot(outlier.alpha = .5, outlier.shape=16, outlier.size=2, outlier.color="grey20", color="grey20", lwd=0.6, fill="grey80", fatten = 1.3)+coord_flip()+
theme_bw()+
ylab("Increased Mean\nSquared Error (%)")+
ggtitle("Total Organic Carbon")+
theme(panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1.1),
axis.title = element_text(size=30),
axis.text = element_text(size=24, colour="black"),
plot.title = element_text(size=30, face = "bold"),
plot.margin = ggplot2::margin(10,30,10,10),
legend.position="none",
axis.title.y = element_blank(),
axis.title.x = element_text(margin = ggplot2::margin(t = 10, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
## Warning: Duplicated aesthetics after name standardisation: size
TOC.Node.plot <- ggplot(TOC.all.rf.import, aes(x = Variable, y = IncNodePurity))+stat_boxplot(geom ='errorbar', width = 0.6, size=0.6, color="grey20", lwd=.1)+ geom_boxplot(outlier.alpha = .5, outlier.shape=16, outlier.size=2, outlier.color="grey20", color="grey20", lwd=0.6, fill="grey80", fatten = 1.3)+coord_flip()+
theme_bw()+
ylab("Increased\nNode Purity")+
ggtitle("Total Organic Carbon")+
theme(panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1.1),
axis.title = element_text(size=30),
axis.text = element_text(size=24, colour="black"),
plot.title = element_text(size=30, face = "bold"),
plot.margin = ggplot2::margin(10,30,10,10),
legend.position="none",
axis.title.y = element_blank(),
axis.title.x = element_text(margin = ggplot2::margin(t = 10, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
## Warning: Duplicated aesthetics after name standardisation: size
Imp.sum <- ggarrange2(Mo.MSE.plot,
Mo.Node.plot,
U.MSE.plot,
U.Node.plot,
prop_eux.MSE.plot,
prop_eux.Node.plot,
TOC.MSE.plot,
TOC.Node.plot,
ncol=2)

ggsave("Figure Sx Variable importance plots for Monte Carlo random forest.pdf", Imp.sum, height=27, width=17)
11. Appendix - Three phase plot evolution for talks
No shading.
Mo.plot.for.sum.no.shade <- ggplot(mo.pdp, aes(x=Age))+
#annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.4,62),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Mo (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
U.plot.for.sum.no.shade <- ggplot(u.pdp, aes(x=Age))+
#annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-.4,62),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("U (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
prop_eux.plot.for.sum.no.shade <- ggplot(prop_eux.pdp, aes(x=Age))+
#annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.01,1.03),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Proportion euxinic")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
TOC.plot.for.sum.no.shade <- ggplot(TOC.pdp, aes(x=Age))+
#annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-0.04,8.3),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("TOC (wt%)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
summary.plot.no.shade <- ggarrange2(Mo.plot.for.sum.no.shade, prop_eux.plot.for.sum.no.shade, U.plot.for.sum.no.shade, TOC.plot.for.sum.no.shade, ncol=2, heights=c(1,1))

ggsave("Figure 2 Partial Dependence Plot no shading.pdf", summary.plot, height=14, width=26)
Neoproterozoic shading.
Mo.plot.for.sum.shade.1 <- ggplot(mo.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.4,62),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Mo (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
U.plot.for.sum.shade.1 <- ggplot(u.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-.4,62),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("U (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
prop_eux.plot.for.sum.shade.1 <- ggplot(prop_eux.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.01,1.03),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Proportion euxinic")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
TOC.plot.for.sum.shade.1 <- ggplot(TOC.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
#annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-0.04,8.3),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("TOC (wt%)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
summary.plot.shade.1<- ggarrange2(Mo.plot.for.sum.shade.1, prop_eux.plot.for.sum.shade.1, U.plot.for.sum.shade.1, TOC.plot.for.sum.shade.1, ncol=2, heights=c(1,1))

ggsave("Figure 2 Partial Dependence Plot with shading part 1.pdf", summary.plot.shade.1, height=14, width=26)
Neoproterozoic and early Paleozoic shading.
Mo.plot.for.sum.shade.2 <- ggplot(mo.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.4,62),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Mo (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
U.plot.for.sum.shade.2 <- ggplot(u.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-.4,62),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("U (ppm)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
prop_eux.plot.for.sum.shade.2 <- ggplot(prop_eux.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_cartesian(xlim=rev(c(298.9,1000)), ylim=c(-.01,1.03),expand=FALSE)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("Proportion euxinic")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top", legend.justification = c(0, 0),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
TOC.plot.for.sum.shade.2 <- ggplot(TOC.pdp, aes(x=Age))+
annotate(geom="rect", xmax=Inf, xmin=541, ymin=-Inf, ymax=Inf, fill="#C9DEE8", alpha=0.2)+
annotate(geom="rect", xmax=541, xmin=418, ymin=-Inf, ymax=Inf, fill="#9DC2D6", alpha=0.4)+
#annotate(geom="rect", xmax=418, xmin=-Inf, ymin=-Inf, ymax=Inf, fill="#619FBB", alpha=0.4)+
geom_ribbon(aes(ymin = min, ymax = max), alpha=1, fill="grey85", color="grey70", size=.3)+
geom_ribbon(aes(ymin = perc.25, ymax = perc.75), alpha=1, fill="grey60", color="grey50", size=.5)+
theme_bw()+
coord_geo(xlim=rev(c(298.9,1002)), expand=FALSE, ylim=c(-0.04,8.3),
pos = as.list(rep("bottom", 1)),
abbrv=list(TRUE),
dat = list(periods.edit),
height = list(unit(2, "lines")),
bord=list(c("left", "bottom", "right")), lwd=as.list(c(1)), size=8)+
scale_x_reverse(breaks=c(300,400,500,600,700,800,900,1000))+
ylab("TOC (wt%)")+xlab("Time (Ma)")+
theme(plot.margin = ggplot2::margin(.1,1,.3,1,"cm"),panel.border = element_rect(fill=NA,color="black", size=2,linetype="solid"),
axis.ticks = element_line(size=1),
axis.line = element_line(lineend = 'square'),
axis.title = element_text(size=34),
axis.text = element_text( size=26, color="black"),
legend.title = element_text(size=20),
legend.text = element_text( size=16),
axis.ticks.length = unit(5, "points"),
legend.position="top",
legend.justification = c(0, 0),
axis.title.y = element_text(vjust=3),
panel.grid.major = element_blank(),panel.grid.minor = element_blank())
summary.plot.shade.2 <- ggarrange2(Mo.plot.for.sum.shade.2, prop_eux.plot.for.sum.shade.2, U.plot.for.sum.shade.2, TOC.plot.for.sum.shade.2, ncol=2, heights=c(1,1))

ggsave("Figure 2 Partial Dependence Plot with shading part 2.pdf", summary.plot.shade.2, height=14, width=26)